Synchrony-optimized Networks of Non-identical Kuramoto Oscillators 
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In this letter we discuss a method for generating synchrony-optimized coupling architectures of 
Kuramoto oscillators with a heterogeneous distribution of native frequencies. The method allows 
us to relate the properties of the coupling network to its synchronizability. These relations were 
previously only established from a linear stability analysis of the identical oscillator case. We 
further demonstrate that the heterogeneity in the oscillator population produces heterogeneity in 
the optimal coupling network as well. Two rules for enhancing the synchronizability of a given 
network by a suitable placement of oscillators are given: (i) native frequencies of adjacent oscillators 
must be anti-correlated and (ii) frequency magnitudes should positively correlate with the degree 
of the node they are placed at. 



An interesting feature of populations of coupled os- 
cillators is synchronization Questions about prop- 
erties of coupling schemes that favor synchronization 
have received much attention in the recent literature 
0, i, i, i, i, S i, i, E El, m • Most of this work has so 
far concentrated on systems of identical oscillators and 
an analysis of the stability of the synchronized state via 
the eigenvalue ratio [l^ and correlations between it and 
properties of the coupling network. Typical findings are 
that network topologies with stable synchronized states 
tend to be very homogeneous in degree- and load dis- 
tributions, very small and characterized by having very 
few short loops and a high degree of disassortative mix- 
ing between node degrees [10|, Indeed, very recent 
studies characterized the most synchronizable networks 
as 'superhomogeneous' [l^. However, such networks are 
rarely found in nature. Part of the homogeneity of the 
networks might be attributable to the fact that all nodes 
are subject to identical dynamics. Variance in the oscilla- 
tor populations is far more likely in nature and seems to 
be the more relevant case in many practical applications. 
Inhomogeneity in the optimal coupling arrangement is 
very likely to arise from the inhomogeneity in the oscil- 
lator population. Thus, understanding the interplay be- 
tween the native frequency of an oscillator and its place 
in the coupling network is a very interesting problem. 

One step in this direction has already been made in 
Ref. p^], where networks of non- identical Kuramoto 
oscillators are rewired by favoring connections between 
oscillators with similar average frequencies. This study, 
however, does not investigate the connection between the 
native frequencies of oscillators and the topology of the 
optimized networks. The latter is the main focus of this 
letter, where some first results about the optimal cou- 
pling of a heterogeneous oscillator population are dis- 
cussed. 

The eigenvalue ratio, which is used in most of the 
previous work, is a measure of the linear stability of 
the synchronized state of coupled identical oscillators. 
Though some network properties that were derived from 
the eigenvalue anlvsis seem to also give rise to better svn- 



chronizability for non- identical oscillators Q , the method 
itsself can make no statement about the non-identical os- 
cillator case. To address the latter issue, we drop the 
generality of the eigenratio analysis and concentrate on 
a specific oscillator dynamics, given by the Kuramoto 
model 
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where the 0i, i — 1, ...,N describe the phases of N oscil- 
lators, which have native frequencies LUi and are coupled 
with strength a via a network described by its adjacency 
matrix A = (ay ). For simplicity we only consider sym- 
metric unweighted coupling. Our choice of the dynami- 
cal model is motivated by the simplicity of the Kuramoto 
model and the fact that it can also be understood as an 
jroximation to a more complex dynamical situation 
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To measure the degree of synchronization of the os- 
cillators, one frequently introduces the order parameter 
r{t) = |l/-/VX]j exp(z0j(i))| Q. For perfect synchro- 
nization, one has r — 1, whereas r ~ 0{N^^/'^) for ran- 
domly drawn omega's and small coupling. A second or- 
der phase transition at some critical coupling strength 
CTc between a desynchronized and a synchronized phase 
is found for mariy coupling schemes [14| . For global cou- 
pling (see e.g. [l5|) and recently for homogeneous and 
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heterogeneous network topologies [If 
synchronization has been well studied. 

The synchronizability of a network is characterized by 
the dependence of r on the coupling strength a. A partial 
ordering might be defined, by calling a network Ai more 
synchronizable than a network A2, if rAi{<y) > M2(c) 
for all a. If rA-i{<Ji) > fA2(o'i) for one cti, but rAj^{cr2) < 
fA2{^2) for another (T2 no statement can be made. 

In what follows, we are interested in the arrangement 
of the couplings that leads to the best synchronizabil- 
ity for a fixed number of links and a given set of native 
frequencies {i^i},fli, where the frequencies oji are drawn 
from a uniform distribution over f— 1, To address 
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this issue we propose a simple optimization scheme. Af- 
ter chosing a value ct*, we start with an Erdos-Renyi 
random graph (ER r.g.) coupling To determine 

a network's synchronizability, we then numerically inte- 
grate Eq.'s (H)) for the initial condition (j>i{t = 0) = and 
calculate the average value of r over AT integration steps 
r{a*) = l/^Tj2j=^f^^ ■r{t,a*), having first allowed for 
Tici timesteps for relaxation. Tj-ci is chosen such that av- 
erage trajectories {r{t > Ti-ci))^^ are stationary for the 
initial network. 

Then, at each iteration of the optimization, wc sug- 
gest a rewiring in the coupling network's configuration. 
By suggesting a move of a randomly chosen link to some 
randomly chosen link vacancy, the connectivity of the 
network is preserved. A new (connected) configuration 
is accepted, if its average fitness r (cr* ) is larger than that 
of the old configuration and rejected otherwise. The pro- 
cedure is repeated till no move was accepted during the 
last 1 = 2 j iterations. 

Three observations justify this somewhat arbitrary 
procedure. First, we find that for a sufficiently large 
choice of cr* a network that has superior synchronizabil- 
ity for homogeneous initial conditions (/)i(t = 0) = 0, i = 
1, TV also has superior average synchronizability if ini- 
tial conditions are drawn from a uniform distribution over 
[— TT, tt] . The reason for this is that the system ([1]) loses 
memory of the initial conditions for large enough cou- 
pling. Second, we note that networks which are better 
synchronizable for one (large enough) choice of a* , also 
prove to be better synchronizable at least for coupling 
strengths larger than a* and generally also for a range 
of coupling strengths below a* . Even more, if for two 
configurations Ai and A2 for which rAi{<J*) > i'A2i'^*)i 
configuration Ai will generally also have an earlier onset 
of synchronization (cf. Fig. [T|). Third, it is found that 
optimal networks obtained for different values of cr* > cr^. 
are very similar. Taking averages over an ensemble of 
initial networks, we indeed find almost identical ^(ct)- 
dependencies for ensembles of networks optimized for 
such different values of cr*. All three observations allow 
to conclude that the networks obtained from the above 
optimization procedure have a superior synchronizabil- 
ity in the sense defined by the above ordering. Contrary 
to the above, a choice of a* much smaller than the crit- 
ical coupling strength for ER r.g.'s only leads to more 
synchronizable networks for the specific choice of initial 
conditions, but does not entail enhanced synchronizabil- 
ity if averaged over an ensemble of initial conditions. The 
reason for this appears to be that memory of the initial 
state is kept during the evolution of the oscillator sys- 
tem and hence the optimization generates networks that 
are specifically adapted to a particular initial distribution 
of phases. Best convergence was generally obtained for 
choices of cr* such that r{a*) « .9 in the initial network. 

We seeded the optimization procedure with different 
initial networks, ranging from scale-free networks con- 
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FIG. 1: Comparison of the dependence r{a) for random and 
several partially optimized networks (the label T gives the 
number of optimization steps performed). Averages are per- 
formed in the following way: first a set of frequencies is uni- 
formly drawn from [—1,1]. Then, for this specific set of 
omega- values, averages over ensembles of initial conditions 
randomly drawn from [— vr, n] are calculated (numerically over 
10 realizations of initial conditions). In a second step, av- 
erages over 100 different sets of native frequencies are per- 
formed. The networks contain N = ICQ nodes and have con- 
nectivity (k) = 3.5. For the optimization we chose a* = .7. 



structed after the Barabasi- Albert model 'i^, ER r.g.'s 
and (almost) regular random graphs. In all cases we find 
very similar optimal configurations which have almost 
identical r(cr)-dependences. Thus, though probably no 
global optimum configuration is reached, the similarity of 
the optima clearly shows how networks change towards 
more synchronizable configurations. 

Although in its aim very similar, our optimization pro- 
cedure differs in some essential points from the one used 
in 12], where rewirings that favour connections between 
elements with similar average frequencies are performed. 
The latter may favor partial synchronization, explaining 
the large level of clustering and the relatively small im- 
provements in r for large coupling strength reported in 

Refs. 
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tion in different networks can be very different. Whereas 
in heterogeneous topologies the fully synchronized state 
tends to grow from a unique core synchronization cluster 
around hubs, the transition to full synchronization in ho- 
mogeneous networks manifests itself by the coalescence 
of many small synchronization clusters of similar size. 
In heterogeneous networks the onset of synchronization 
arises for smaller coupling, but this is traded-off against 
a lesser stability of the fully synchronized state and a 
smaller degree of synchronization r for larger couplings. 
By the choice of a large cr* in the balance between the on- 
set of synchronization and stability of the synchronized 
state more emphasis is put on the latter. 

Below we proceed by investigating correlations be- 
tween network configurations and native frequency ar- 



3 




0.4 0.8 1.2 1.6 



FIG. 2: Structure of an optimized network (a) and comparison 
of the sigma-dependence of the degree of synchronization r be- 
tween the original (open circles) and the optimized network 
(filled circles) (b). In the network illustration, nodes associ- 
ated with negative native frequencies are drawn in red, others 
in green. Connections between negative-positive omega pairs 
are drawn as solid lines, others as fat dashed lines. Note, that 
there are only two cases where nodes with frequencies of the 
same sign are neighbours. 
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rangements generated by the optimization approach. To 
gain intuition, Fig. [^h- displays a small optimal network 
and the corresponding assignment of native frequencies. 
Closer inspection of the optimal network strongly sug- 
gests that its structure is characterized by a pairing of 
the native frequencies, such that nodes with uj > are 
always surrounded by other nodes with w < and vice 
versa. We introduce two quantities to describe these cor- 
relations between the w-values of adjacent nodes. First, 
P- measures the proportion of links connecting nodes 
with omegas of different signs. Second, a correlation co- 
efficient 

where (uj) = ^/N'^^uji additionally quantifies how often 
nodes with large negative omega are adjacent to nodes 
with strongly positive omega. For arbitrary assignments 
of native frequencies one has — 1/2 and (c^j) — 0. 
We proceed by investigating the dependence of p_ and 
as well as of several network characteristics on syn- 
chronizability. For this, we determine averages of the 
clustering coeficient, pathlength, diameter, assortative- 
ness, load 

(i) 

where d ^ ^ = 1 if the shortest path from node j to node 

k passes through node i and d^- ^ =0 otherwise, load 
and degree variance and maximum load during the net- 
work optimization. In the following, we limit the dis- 
cussion to changes in the network, where the starting 
point of the optimization was a random graph. Start- 
ing with other initial conditions, e.g., scale-free networks, 



FIG. 3: Dependence of the clustering coefficient (a), den- 
sity of 4- loops (b), assortativeness (c), maximum load (d), 
load variance (e), average pathlength (f), degree variation (g), 
fraction of plus-minus pairs (h) and of the correlation of ad- 
jacent native frequencies (i) on synchronizablity measured by 
r((j*) during the optimization. The data are for a network of 
100 nodes with connectivity (fc) = 3.5. The initial network 
corresponded to an ER r.g., for other initial conditions like 
scale-free networks or regular random graphs similar trends 
are found. 



changes in most network properties stem from a reorga- 
nization of the hub structures towards more 'democratic' 
arrangements. However, in the last stages of the evolu- 
tion, trends are very similar to those we discuss for ER 
r.g.'s below. 

Our most important finding is the observation that 
the optimal networks are characterized by a strong anti- 
correlation of adjacent positive and negative frequencies. 
Both, the number of plus-minus pairs p- as well as the 
anti-correlation are found to increase strongly during the 
evolution of the networks towards better synchronizabil- 
ity (cf. Fig. [3)i,i). Further, on the one hand, the changes 
in the assortativeness, cliquishness, load, maximum load 
and load variance as well as in the average pathlengths 
and diameter (not shown) corroborate results obtained 
from the eigenratio analysis for identical oscillators. Het- 
erogeneity in the oscillator population does not appear to 
affect network properties that are related to the transmis- 
sion of information. On the other hand, the network con- 
figurations with optimal synchroniz ability are not found 
to be the most homogeneous ones, and while exhibit- 
ing fewer triangles, are not generally marked by lower 
densities of short loops than random networks. Since 
triangles always involve at least one adjacent pair of os- 
cillators with native frequencies of the same sign, the 
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decay in the triangle density, and indeed the density of 
odd loops generally, during optimization follows from the 
increase in the number of frequency-anticorrelated node 
pairs. Hence the triangle density can not be used as an 
indicator for general short loop densities. For this rea- 
son, we also measured the density of loops of length four, 
which is indeed found to increase during the optimization 
(cf. Fig.lSb). 

Next, we discuss the development in the networks' de- 
gree variance during their evolutions. For this, it appears 
worthwhile to observe that the largest degree of anti- 
correlation between adjacent native frequencies can be 
realized in very inhomogeneous networks, where frequen- 
cies of very large (or very small) magnitude attract more 
links than frequencies of average magnitude. This lets 
one understand that degree homogeneity and enhanced 
anti-correlation of adjacent node frequencies are conflict- 
ing demands. 

Figure[3|; displays a minimum in the dependence of the 
degree variance on synchronizability. Initially, the degree 
variance is reduced, but around r = 0.98 this trend is re- 
versed. We hypothesize that this point is associated with 
an arrangement of the network where no further increase 
in the strength of anti-correlations between adjacent os- 
cillators can be found in very regular coupling schemes. 
Instead, the degree of anticorrelation can only be further 
enhanced by placing native frequencies of large absolute 
value on nodes of higher degrees, which gives them ac- 
cess to many neighbours of low frequency magnitude. As 
oscillators synchronize with the mean frequency {lo) w 
(for even distributions of the w's), this placement can be 
understood as resulting in a strong drag on large frequen- 
cies towards the mean frequency. The same effect in the 
correlations could be achieved by placing frequencies of 
low magnitude on the higher degree nodes and surround 
them by frequencies of large magnitude. However, the 
latter arrangement clearly favours the build-up of small 
clusters of oscillators synchronizing with frequencies of 
large magnitude. 

In fact, the optimal network displayed in Fig[T^ ap- 
pears very close to the first arrangement of the native 
frequencies. For example, both frequencies of the largest 
magnitude w = — .92 and lu = —.91 are assigned to the 
two nodes with the most connections, while the next two 
belong to nodes with the second most connections, etc. 
We also measured the dependence of the average magni- 
tude of the native frequency on the node degree (|Li;|)(/c) 
and find a steady increase from (|a;|)(l) w .15 for nodes 
of degree one to (|ci;|)(8) « .75 for the highest degree 
node of the optimal network. Corresponding with this, 
the average magnitude of the difference between the fre- 
quency of a node and the average frequency of its near- 
est neighbours Awn.n. = \w — l/^i X^i '^y'^il grows from 
(Awn.n.)(l) ~ .7 to (Au;„.„.)(8) « 1.08. Clearly, the opti- 
mal network configuration for heterogeneous oscillators is 
not 'superhomogeneous'. Instead 'hub' nodes are found 
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FIG. 4: Comparison of the synchronizability of (a) ER and 
(b) scale-free networks with a random distribution of native 
frequencies (open squares) and distributions with strong anti- 
correlations between adjacent frequencies (filled circles) and 
frequency arrangements with an additional correlation be- 
tween degree and frequency magnitude (open circles). The 
data are for networks of size = 1000 with (k) = 6. 



which are distinguished by native frequencies of large 
magnitude, whereas low degree nodes are typically as- 
sociated with frequencies of relatively small magnitude. 

To obtain another independent confirmation for the 
main contention of the paper, that anti-correlated ar- 
rangements of native frequencies on nodes strongly en- 
hance the synchronization behaviour, we proceed by com- 
paring the r(cr)-dependence of networks with random and 
strongly anti-correlated w-assignments. The latter can 
easily be generated by proposing frequency swaps be- 
tween nodes on a fixed network and accepting them if 
they lead to a decrease of c^j. Figure [4] shows the de- 
pendence of r on (7 for scale-free and ER networks. The 
data corroborate our assertion. It is also quite interest- 
ing to generate frequency distributions that additionally 
positively correlate frequency magnitude and node de- 
gree. In the same way as above such an architecture 
can be created by conditionally swapping frequencies on 
nodes and minimizing for Ccj — c\^\k, where c\t^\k is a stan- 
dard correlation coefficient between the \uj\ values and 
the node degrees k. Figure [4] shows that this frequency 
arrangement leads to a lower order parameter r for low 
coupling, but increased r for high cr's. This behaviour 
is straightforward to understand from our previous ar- 
guments. The additional correlation of and the node 
degree entails larger differences between native frequen- 
cies belonging to a hub and to its average neighbour. As 
hubs typically serve as germs for clusters of synchrony 
IgL 17 1, the latter are harder to form in this scenario. 

In summary, we have discussed a method for generat- 
ing networks of non-identical Kuramoto oscillators with 
enhanced synchronizability. Using this method, some 
correlations between properties of the coupling architec- 
ture and synchronizability which were previously only 
derived for the identical oscillator case could be demon- 
strated for an example system of non-identical oscillators. 
As another important point, evidence has been given, 
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that the synchronizability depends not only on the topol- 
ogy of the coupling architecture, but is also strongly in- 
fluenced by the assigncmcnt of frc'qucncics to nodes. Wo 
have demonstrated that the synchronization behaviour 
can be strongly enhanced by a frequency arrangement 
that is marked by strong anti-correlations between fre- 
quencies of adjacent nodes. A trade-ofi^ between this ef- 
fect and the requirements for good communication leads 
to some degree of heterogeneity in the networks with opti- 
mal synchronization properties. Further, we have shown 
that the optimal assignment is such that nodes with more 
than an average number of connections are assigned na- 
tive frequencies of large magnitude, while frequencies of 
small magnitude typically belong to small degree nodes. 
Our findings give a set of rules for constructing more syn- 
chronizable networks of heterogeneous oscillators. These 
might be of practical importance for engineering prob- 
lems, but also give predictions that should be tested for a 
range of biological networks where synchronization plays 
a role. 
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